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ABSTRACT 

This report describes the physical models employed in GPSOMC, the modeling module of the 
GIPSY software system developed at JPL for analysis of geodetic Global Positioning Satellite (GPS) 
measurements. Det ails of the various contributions to range and phase observables are given, as well 
as the partial derivatives of the observed quantities with respect to model parameters. A glossary 
of parameters is provided to enable persons doing data analysis to identify quantities in the current 
report with their counterparts in the computer programs. The present version is the second revision 
of the original document, JPL Publication 87-21, dated September 15, 1987, which it supersedes. The 
modeling is expanded to provide the option of using Cartesian station coordinates; parameters for the 
time rates of change of universal time and polar motion are also introduced. 
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SECTION 1 


INTRODUCTION 

In interpreting measurements of range from satellites of the Global Positioning System (GPS) 
to ground-based receivers, the observables are passed through a multiparameter estimation routine 
(“filter*) to estimate the parameters of a model. This model describes the spacecraft orbits and the 
motions of the Earth-fixed receivers and supplies to the filter a priori values of computed observables, 
as well as their partial derivatives with respect to model parameters. During 1984-85, software 
was developed at JPL to perform these functions. The purpose of the present report is to describe 
that portion of the software which concerns the modeling of receiver locations, motions of the whole 
Earth, and computation of observables and their partial derivatives. Modeling of satellite orbits 
and parameter estimation form separate units in the software chain and are described in separate 
documents. 

Many aspects of the model required to describe GPS range measurements are identical to the 
modeling developed during the past decade for Very Long Baseline Interferometry (VLBI). Conse- 
quently much of the content of this document, as well as of the associated software package GPSOMC, 
is borrowed from the VLBI modeling and parameter estimation package MASTERFIT and the doc- 
ument giving the details of its models (Sovers and Fanselow, 1987). There are two major differences. 
First, in satellite range measurements, the sources (transmitters) are not at infinity, and they contain 
time standards of their own. Together with the need for orbit determination, this makes the situation 
considerably more complicated than the VLBI case, where signals are received from fixed sources 
at infinite distances. The second, minor, model difference is the lack of necessity of describing the 
passage of the signal through the Earth’s ionosphere. It is assumed that ionospheric effects will always 
be removed by performing measurements at two frequencies. 

Three major model components which will be discussed here are geometry, clocks, and tropo- 
sphere. Section 2 deals with the coordinate frame for the model and establishes methods for calcu- 
lating the position of the receiver in that frame, employing the best current models of whole Earth 
motions and local tidal deformations. This section is nearly identical to the corresponding coordi- 
nate frame description in the MASTERFIT document, with slight differences reflecting differences 
in the implementation of minor effects. The next section (Sec. 3) defines the observables and the 
intimately associated models of behavior of space vehicle and receiver clocks. Section 4 presents the 
model employed to describe the passage of the signal through the troposphere. All the partial deriva- 
tives of observables with respect to model parameters are given in Sec. 5. Values of the physical 
constants used in the GPSOMC software, as well as a complete list of the parameters presently available 
for adjustment, are given in the appendices. 

As of June 1990, the model description in this report corresponds to the Fortran code used to 
generate the executable GPSOMC . EXE ; 128 on the Section 335 VAX computers. This correspondence 
holds for the dominant components of the observable models. Some of the code which implements 
aspects of modeling that are seldom used ( e.g. } UTPM rates, nutation parameters, time-dependent 
station positions and zenith tropospheres) has not yet been thoroughly checked; it is intended that 
complete testing will be performed in the near future. 
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SECTION 2 


COORDINATE FRAME AND GEOMETRY 

The geometric range is that portion of the distance between a satellite transmitter and an Earth- 
fixed receiver which would be measured by perfect instrumentation, and perfectly synchronized, if there 
were a perfect vacuum between the transmitter and receiver. It is by far the largest component of the 
observed range. The main complexity of this portion of the model arises from the numerous coordinate 
transformations necessary to relate the inertial reference frame used for locating the spacecraft to the 
Earth-fixed reference frame in which station locations are represented. 

In the following we will assume that “geocentric inertial reference frame” means a geocentric, 
equatorial frame with the equator and equinox of J2000 as defined by the 1976 IAU (International 
Astronomical Union) conventions with the 1980 nutation series (Seidelmann, 1982, and Kaplan, 1981). 
On the other hand, we will mean by “Earth-fixed reference frame” some reference frame tied to 
the mean surface features of the Earth. This is a right-handed version of the CIO (Conventional 
International Origin) reference system with the pole defined by the 1903.0 pole. Implementation is 
accomplished by defining the position of one of the fiducial stations, and measuring the positions 
of the other receivers. This section provides the details for the evaluation of the geometric range, 
including receiver coordinates, tidal effects, and the transformation from Earth-fixed to geocentric 
inertial coordinate systems. 

2.1 RECEIVER LOCATION TIME DEPENDENCE 

Normally the receiver position vector te 0 and its components in the Earth- fixed reference frame 
r 9pi A, z (radius off the spin axis, longitude, and height above the equator, respectively) are considered 
to be time- invariant. An alternative formulation introduces the time rates of change of the station 
coordinates as adjustable parameters. The model is linear, with the components of r£ 0 (t) at time t 
expressed as 

r »p = r % + - to) (2.1) 

A = A 0 + A(t - t 0 ) (2.2) 

z = z° + z(t - t 0 ) (2.3) 

Here t 0 is a reference epoch, at which the receiver coordinates are (r° p , A 0 , z°). Alternatively, a 
Cartesian Earth-fixed frame may be used, with components of r£ 0 (f) expressed as 

X = Xo + X(t — to) (2.4) 

Y = Y 0 + Y(t- t 0 ) (2.5) 

Z = Zo 4- Z(t — to) (2.6) 

The receiver “position vector” may include an antenna phase center offset and an offset to a standard 
benchmark or monument (see Sec. 2.9). In the model development that follows, r Eo includes these 
offsets, but is referred to as the “receiver vector.” Besides the linear time variation of Eqs. (2.4—2. 6), 
modeling allows for estimation of a new benchmark position for every observing session. 

2.2 TIDAL EFFECTS 

In calculating the geometric range, we need to consider the effects of crustal motions on receiver 
locations. Among these deformations are solid Earth tides, tectonic motions, and alterations of the 
Earth’s surface due to local geological, hydrological, and atmospheric processes. In the standard 
Earth- fixed coordinate system, tidal effects modify the receiver location Te 0 by an amount 

^ ~ ^*ol Aocn + Apai (2*7) 

where the three terms are due to solid Earth tides, ocean loading, and pole tide, respectively. Other 
Earth-fixed effects would be incorporated by augmenting the definition of A. 
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2.2.1 Solid Earth Tides 


Alteration of the positions of the receivers by solid Earth tides is rather complicated due to their 
coupling with the ocean tides and to the effects of local geology. We gloss over these complications, and 
employ the simple quadrupole response model described by J. G. Williams (1970), who used Melchior 
(1966) as a reference. Let R # be the position vector of a perturbing source in the Earth-fixed reference 
system, and let te 0 be the receiver position vector in the same coordinates. If h and l are the Love 
numbers, rp a phase shift of the tidal effects, and r p the phase-shifted receiver location vector in the 
Earth- fixed reference system, then the vector of tidal displacements in a local Cartesian frame (x axis 
vertical, y axis eastward, and % axis northward on a spherical Earth) is 

* = to«> to#| r ( 2 - 8 ) 


where 


_ 3fi,r p [ ( r p R ») 2 _ r p R > 

9u R* [ 2 6 

92, — ~ (r p • R,)(y,x p — X,y p )—j 




( 2 . 10 ) 


93. = (r p • R.) y/x l + y* Z, — p = {x p X, + y p Y t ) 
* L y x p + y p 


( 2 . 11 ) 


Here fi a is the ratio of the mass of the disturbing object, s, to the mass of the Earth, and 

R. = [X.,Y„Z.] t (2.12) 

is the vector from the center of the Earth to that body. The summation is over all tide-producing 
bodies, of which we include only the Sun and the Moon. 

The phase-shifted receiver vector is calculated employing a phase lag, or equivalently, employing 
a right-handed rotation, L, through am angle rp about the Z axis of date, r p = Lr^ 0 . This lag matrix, 
L, is: 

( cos rp sin rp 0 A 

— sin rp cos rp 0 j (2.13) 

0 0 l) 

By a positive value of rp we mean that the peak response on an Earth meridian occurs at a time 
St = rp/wE after that meridian containing r^ 0 crosses the tide-producing object, where we is the 
angular rotation rate of the Earth. In the vertical component, the peak response occurs when the 
meridian containing r p also includes R,. The difference between geodetic and geocentric latitude can 
affect this model on the order of the (tidal effect) /(flattening factor) « 0.1 cm. 

To convert the locally referenced strain, £, to the Earth-fixed frame, two additional rotations 
must be performed. The first, W , rotates the vector by an angle, <p Si about the y axis to an equatorial 
system. The second, V , rotates about the resultant z axis by angle, —A,, to bring the displacements 
into the standard Earth-fixed coordinate system. The result is 


where 


A.oj = VW6 

(2.14) 

cos <p 9 0 — sin <P a \ 

0 10 
sin <p a 0 cos <p a J 

(2.15) 
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and 


V = 


C08 A 
8 in A, 
0 


(2.16) 


— sin A, O' 
cos Aj, 


Actually, the product of these two matrices is coded: 

( cos A, cos <f> 9 — sinA # —cos A 

sin A, cos <j> 9 cos A, —sin A 
sin <t> 9 0 cos 


« sin <f> 9 \ 
9 sin <j> 9 

» 4 > - J 


(2.17) 


2.2.2 Ocean Loading 

This section is concerned with one of the secondary tidal effects, t.e., the elastic response of the 
Earth’s crust to ocean tides, which move the receivers to the extent of a few cm. Such effects are 
commonly labeled “ocean loading.” The most complete recent model appears to be that described 
by Pagiatakis (1982, 1988) and by Pagiatakis, Langley, and Vanicek (1982), which is implemented 
as an option in the current GPSOMC. It was recently generalized to include effects of anisotropy and 
viscoelasticity of the Earth; differences of displacements with those given by the 1982 model are well 
below 1 cm. The formulation and coding are general enough, however, to permit other inputs to 
be used in place of the Pagiatakis ocean loading model. Because the receiver motions caused by 
response to ocean tides appear to be limited to approximately 3 cm for sites well removed from the 
coast, no estimation capability was deemed necessary at present. This decision is supported by the 
fact that for locations near the coast, where the effects may be more sizable, and which would thus 
be expected to produce data useful in parameter estimation, the elastic response modeling is as yet 
inadequate (Agnew, 1982). The present model entails deriving an expression for the locally referenced 
displacement 6 due to ocean loading. In the vertical, N-S, E-W local coordinate system at time t, 

N 

Sj co 8 ( W< t + Vi - S’) (2.18) 

i= 1 

The quantities w* (frequency of tidal constituent t) and V 9 (astronomical argument of constituent t) 
depend only on the ephemeris information (positions of the Sun and Moon). On the other hand the 
amplitude £? and Greenwich phase lag 6* of each tidal component j are determined by the particular 
model assumed for the deformation of the Earth. As of November 1988, software for calculating 
these deformations at an arbitrary point on the Earth’s surface exists at JPL only for the Pagiatakis- 
Langley model. Six tidal components are included [N = 6 in Eq. (2.18)]: the S 2 , Hi, O*, N 2 , 
and Pi tides, all of which have periods close to either 12 or 24 hours. The local displacement vector is 
transformed via Eqs. (2.17) and (2.14) to the displacement A ocn in the standard Earth-fixed frame. 

Input to GPSOMC provides for specification of an arbitrary number of frequencies and astronomical 
arguments w* and V*, followed by tables of the local distortions and their phases, £? and 6* 9 calculated 
from the ocean tidal loading model of choice. In particular, longer-period tidal constituents can be 
accommodated in this fashion. 

There are presently four choices of models. As mentioned above, the six components of the 
Pagiatakis- Langley model can be calculated by separate software for an arbitrarily located receiver 
to generate input tables for GPSOMC. This has been done for all stations commonly employed in JPL 
VLBI experiments. There is no comparable facility to obtain amplitudes and phases for the Agnew 
(1982), Scherneck (1983), or Goad (1983) models. Consequently for these three models GPSOMC input 
tables only exist for the limited set of stations considered by these authors. Agnew considers only five 
components (omitting Pi), while Scherneck includes five components in addition to the Pagiatakis- 
Langley set: Q 1 , A//, M mt and S 9a . These all have amplitudes of 1 mm or smaller and, thus, 

are not expected to be significant at the present level of experimental accuracy. Goad’s MERIT 
standard model only specifies vertical displacements, but includes three components (K 2 , Q 1 and 
M/) in addition to Pagiatakis-Langley. Due to their bulk, none of the tables of tidal amplitudes are 
reproduced here, but are available on request in computer-readable form. 
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2.2.3 Pole Tide 


In addition to the effects of ocean tides, another secondary tidal effect is the displacement of a 
receiver by the elastic response of the Earth’s crust to shifts in the spin axis orientation. The spin axis 
is known to describe a circle of « 20-m diameter at the north pole. Depending on where the spin axis 
pierces the crust at the instant of a measurement, the “pole tide* displacement will differ from time 
to time. This effect must be included if centimeter accuracy is desired, especially for measurements 
spanning an appreciable fraction of a year. 

Yoder (1984) derived an expression for the displacement of a point at latitude longitude A due 
to the pole tide: 


6 


w\R 


[sin <f> cos 4>(x cos A + y sin A) hr 


+ cos 2<f>(x cos A + y sin A) 

+ sin ^(—x sin A + ycos A)ZA] 


(2.19) 


Here we is the rotation rate of the Earth, R the radius of the (spherical) Earth, g the acceleration 
due to gravity at the Earth’s surface, and h and l the customary Love numbers. Displacements of the 
spin axis from the 1903.0 CIO pole position along the x and y axes are given by x and y. Eq. (2.19) 
shows how these map into receiver displacements along the unit vectors in the radial (r), latitude (^), 
and longitude (A) directions. With the standard values we — 7.292 X 10“ 6 rad/sec, R — 6378 km, 
and g = 980.665 cm/sec 2 , the factor a>| R/g = 3.459 X 10“ 3 . Since the maximum values of x and y 
are of the order of 10 meters, and h « 0.6, l » 0.08, the maximum displacement due to the pole tide 
is 1 to 2 cm, depending on the location of the receiver (^, A). 

As in the case of the previously considered tidal effects, the locally referenced displacement 6 is 
transformed via the transformation Eq. (2.17) to give the displacement A po j in the standard Earth- 
fixed coordinate system. After each of the locally referenced tidal displacements has been transformed 
to these coordinates, the receiver location is 


te = r^ 0 + + Aocn + A p oi 


( 2 . 20 ) 


2.3 TRANSFORMATION FROM EARTH-FIXED TO GEOCENTRIC 
INERTIAL COORDINATE SYSTEMS 

The Earth is approximately an oblate spheroid, spinning in the presence of two massive moving 
objects (the Sun and the Moon) which are positioned such that their time-varying gravitational effects 
not only produce tides on the Earth but also subject it to torques. In addition, the Earth is covered 
by a complicated fluid layer, and also is not perfectly solid internally. As a result, the orientation 
of the Earth is a very complicated function of time, which to first order can be represented as the 
composite of a time-varying rotation rate, a wobble, a nutation, and a precession. The exchange of 
angular momentum between the solid Earth and the fluids on its surface is not readily predictable, 
and thus must be continually determined experimentally. Nutation and precession are well modeled 
theoretically. At the centimeter level, however, even these models are not completely adequate. 

The rotational transformation, Q , of coordinate frames from the Earth- fixed frame to the geocen- 
tric inertial frame is composed of 6 separate rotations (actually 10, since the nutation transformation, 
N t consists of 3 transformations in itself, as does the “perturbation* transformation, O) applied to a 
vector in the Earth- fixed system: 

Q = QPNUXY (2.21) 

Thus, if te is a receiver location expressed in the Earth-fixed system, e.y., the result of Eq. (2.20), 
that location, r /, expressed in the geocentric inertial (J2000) system is 

rj = Qt e (2.22) 
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Note that since we rotate the Earth rather than the celestial sphere, our rotation matrices, fi, P, and 
N , will be the transposes of those used to rotate the inertial system of J2000 to an inertial system of 
date. 

2.4 UT1 AND POLAR MOTION 


The first transformation, F, is a right-handed rotation about the x axis of the Earth-fixed frame 
by an angle 63. Currently, the Earth-fixed frame is the 1903.0 CIO frame, except that the positive 
y axis is at 90 degrees east (Moscow). The x axis is coincident with the 1903.0 meridian of Greenwich, 
and the % axis is the 1903.0 standard pole. The Y rotation matrix is 

(10 0 A 

Y — 10 cos ©3 sin ©3 I (2.23) 

Y 0 — sin ©3 cos ©2 J 

where ©2 is identical to the © y that may be readily obtained from the pole position published by 
the International Earth Rotation Service (IERS, 1989) or the International Radio Interferometric 
Surveying (IRIS) project [I AG (International Association of Geodesy), 1987], 

The next rotation in sequence is the right-hand rotation through an angle ©1 about the y axis 
obtained after the previous rotation has been applied: 



(2.24) 


In this rotation, ©1 is identical to the ©* obtainable from the IERS or IRIS value for the pole 
position. Note that we have incorporated in the matrix definitions the transformation from the left- 
handed system used by IERS to the right-handed system we use. Note also that any source of polar 
motion data can be used provided it is represented in a left-handed system. The only effect would be 
a change in the definition of the Earth-fixed reference system. 

The application of m XY* to a vector in the Earth-fixed system of coordinates expresses that 
vector as it would be observed in a coordinate frame whose 2 axis was along the Earth’s ephemeris 
pole. The third rotation, U, is about the resultant s axis obtained by applying U XY.” It is a rotation 
through the angle, -H } where H is the hour angle of the true equinox of date (t.e., the dihedral angle 
measured westward between the xa plane defined above and the meridian plane containing the true 
equinox of date). The equinox of date is the point defined on the celestial equator by the intersection 
of the mean ecliptic with that equator. It is that intersection where the mean ecliptic rises from below 
the equator to above it (ascending node). 


( cos H — sin H 0\ 
sin if cos if 0 I 

0 0 l) 


(2.25) 


The angle H is composed of two parts: 

H = h 1 + a E (2.26) 

where is the hour angle of the mean equinox of date, and a E is the difference in hour angle of 
the true equinox of date and the mean equinox of date, a difference which is due to the nutation of 
the Earth. This set of definitions is cumbersome and couples the nutation and precession effects into 
Earth rotation measurements. However, in order to provide a direct estimate of conventional UT 1 
(universal time) it is convenient to endure this historical approach, at least for the near future. 

As of 1989, the software provides the option of estimating UTPM rates. The defining expressions 
for rates of x, y polar motion and UT 1 are 

©1 = ©1 + ©i(* — *re/) 

©2 = ©2 + © 2 (* — *re/) 

H x = Ho + H(t - t ref ) 
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(2.27) 

(2.28) 
(2.29) 



UT 1 is defined to be such that the hour angle of the mean equinox of date is given by the following 
expression (Aoki et &!., 1982, and Kaplan, 1981): 


where 


h . 7 = UT 1 + e h 41 m 50*. 54841 + 8640184V8 12866 T u 
+ 0' .093 104 Tl - 6'. 2 x 10 -6 7^ 

( Julian UTl date) - 2451545.0 


r« = 


36525 


The actual equivalent expression which is coded is 


h 7 =2 jt(C/T1 Julian day fraction) + 673 10*. 54841 

+ 8640184*. 8 12866 T u + 0* .093104 2? - 6* .2 x 10 -6 T* 


(2.30) 

(2.31) 

(2.32) 


This expression produces a time, UT 1 , which tracks the Greenwich hour angle of the real Sun to 
within 16 m . However, it really is sidereal time, modified to fit our intuitive desire to have the Sun 
directly overhead at noon on the Greenwich meridian* Historically, differences of UT 1 from a uniform 
measure of time, such as atomic time, have been used in specifying the orientation of the Earth. Note 
that this definition has buried in it the precession constant since it refers to the mean equinox of date. 

By the very definition of “mean of date” and “true of date,” nutation causes a difference in the 
hour angles of the mean equinox of date and the true equinox of date. This difference, called the 
“equation of equinoxes,” is denoted by ct£ and is obtained accordingly: 


where the vector 



(2.33) 


(2.34) 


is the unit vector, in true equatorial coordinates of date, toward the mean equinox of date. In mean 
equatorial coordinates of date, this same unit vector is just (1,0, 0) T . The matrix is just the 

inverse (or equally, the transpose) of the transformation matrix N which will be defined below to 
effect the transformation from true equatorial coordinates of date to mean equatorial coordinates of 
date. 

Depending on the smoothing used to produce the a priori UTl — UTC series, the short-period 
(t < 35 days) fluctuations in UTl due to changes in the latitude and size of the mean tidal bulge may 
or may not be smoothed out. Since we want as accurate an a priori as possible, it may be necessary 
to add this effect to the UTl a priori obtained from the series, UTl smoot hed • If this option is selected, 
then the desired a priori UTl is given by 


UTl a priori — UTl $rnoo thcd AUTl (2.35) 

UTl Bmoo thed represents an appropriately smoothed a priori measurement of the orientation of the 
Earth (e.#., IERS Bulletin A smoothed or, even better, UTlR) t for which the short period (t < 35 
days) tidal effects either have been averaged to zero, or, as in the case of UT1R } removed before 
smoothing. This AUT1 can be represented as 


A UT1 = 



(2.36) 


where N is chosen to include all terms with a period less than 35 days. There are no other con- 
tributions until a period of 90 days is reached. However, these long-period terms are included by 
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the measurements of the current Earth-orientation measurement services. The values for kij and Ai , 
along with the period involved are given in Table I. The a * for t = 1, 5 are just the angles defined 
below in the nutation series a s /, i', F, D, and ft respectively. 

It is convenient to apply tt UXY* as a group. To parts in 10 12 , XY — YX. However, with the 
same accuracy UXY ^ XYU. Neglecting terms of O(0 2 ) (which produce receiver location errors of 
approximately 6 X 10“ 6 meters): 

( cos H — sin H 
sin H cos H 
sin 0i —sin 02 

Table I 

Periodic Tidally Induced Variations in UTl 
With Periods Less Than 35 Days 


— sin 0i cos fT — sin 02 sin H \ 

— sin ©i sin H -f sin 02 cos H j (2.37) 


Index 

Period 

( day *) 

Argument coefficient 
kii ki2 ki 3 ki4 fcjft 

Ai 

( O '. OOOl ) 

1 

5.64 

1 

0 

2 

2 

2 

- 0.02 

2 

6.85 

2 

0 

2 

0 

1 

- 0.04 

3 

6.86 

2 

0 

2 

0 

2 

- 0.10 

4 

7.09 

0 

0 

2 

2 

1 

- 0.05 

5 

7.10 

0 

0 

2 

2 

2 

- 0.12 

6 

9.11 

1 

0 

2 

0 

0 

- 0.04 

7 

9.12 

1 

0 

2 

0 

1 

- 0.41 

8 

9.13 

1 

0 

2 

0 

2 

- 0.99 

9 

9.18 

3 

0 

0 

0 

0 

- 0.02 

10 

9.54 

-1 

0 

2 

2 

1 

- 0.08 

11 

9.56 

-1 

0 

2 

2 

2 

- 0.20 

12 

9.61 

1 

0 

0 

2 

0 

- 0.08 

13 

12.81 

2 

0 

2 

-2 

2 

0.02 

14 

13.17 

0 

1 

2 

0 

2 

0.03 

15 

13.61 

0 

0 

2 

0 

0 

- 0.30 

16 

13.63 

0 

0 

2 

0 

1 

- 3.21 

17 

13.66 

0 

0 

2 

0 

2 

- 7.76 

18 

13.75 

2 

0 

0 

0 

-1 

0.02 

19 

13.78 

2 

0 

0 

0 

0 

- 0.34 

20 

13.81 

2 

0 

0 

0 

1 

0.02 

21 

14.19 

0 

-1 

2 

0 

2 

- 0.02 

22 

14.73 

0 

0 

0 

2 

-1 

0.05 

23 

14.77 

0 

0 

0 

2 

0 

- 0.73 

24 

14.80 

0 

0 

0 

2 

1 

- 0.05 

25 

15.39 

0 

-1 

0 

2 

0 

- 0.05 

26 

23.86 

1 

0 

2 

-2 

1 

0.05 

27 

23.94 

1 

0 

2 

-2 

2 

0.10 

28 

25.62 

1 

1 

0 

0 

0 

0.04 

29 

26.88 

-1 

0 

2 

0 

0 

0.05 

30 

26.98 

-1 

0 

2 

0 

1 

0.18 

31 

27.09 

-1 

0 

2 

0 

2 

0.44 

32 

27.44 

1 

0 

0 

0 

-1 

0.53 

33 

27.56 

1 

0 

0 

0 

0 

- 8.26 

34 

27.67 

1 

0 

0 

0 

1 

0.54 

35 

29.53 

0 

0 

0 

1 

0 

0.05 

36 

29.80 

1 

-1 

0 

0 

0 

- 0.06 

37 

31.66 

-1 

0 

0 

2 

-1 

0.12 

38 

31.81 

-1 

0 

0 

2 

0 

- 1.82 

39 

31.96 

-1 

0 

0 

2 

1 

0.13 

40 

32.61 

1 

0 

-2 

2 

-1 

0.02 

41 

34.85 

-1 

-1 

0 

2 

0 

- 0.09 
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2.5 NUTATION 


With the completion of the UT 1 and polar motion transformations, we are left with a receiver 
location vector, Ydate ■ This is the receiver location in true equatorial inertial coordinates of date. The 
last set of transformations are nutation, AT, precession, P, and the perturbation rotation, ft, applied in 
that order. These transformations give the receiver location, r /, in geocentric inertial J2000 equatorial 
coordinates: 

r/ = nPNTdau (2.38) 

Both the nutation and precession rotation angles are defined relative to their values at Julian date 
2451545.0 (J2000). The angles are computed from trigonometric and polynomial series as a function 
of Barycentric Dynamic Time (TDB, Temps Dynamique Barycentrique). This time scale, which is 
also used to reference the planetary ephemeris, is related to Terrestrial Dynamic Time (TDT, Temps 
Dynamique Terrestre) by (Kaplan, 1981) 

TDB = TDT + 0 J .001658 sin (g + 0.0167sin^) (2.39) 

g = 2*(357°.578 + 35999°.050 TDT)/ 360° (2.40) 

The time scale TDT runs at the same rate as International Atomic Time ( TAI , Temps Atomique 
International) and is offset from TAI by a defined constant, 

TDT = TAI + 32M84 (2.41) 


All time arguments used in nutation and precession computations are measured in centuries of TDB 
from J2000 [c/. Eq. (2.31)]. 

The transformation matrix N is a composite of three separate rotations (Melbourne et al., 1968): 

1. A(e) : true equatorial coordinates of date to ecliptic coordinates of date. 

(10 0 \ 

A(e) = I 0 cos c sine I (2.42) 

y 0 — sin e cos e J 

where e is the true obliquity of the ecliptic. 

2. C 7 (£V>) : nutation in longitude from ecliptic coordinates of date to mean ecliptic coordinates of 
date. 


( cos Sip sin Sip (A 

— sin Sip cos Sip 0 I (2-43) 

0 0 1 J 

where Sip is the nutation in ecliptic longitude. 

3. A r (e) : ecliptic coordinates of date to mean equatorial coordinates. 

In ecliptic coordinates of date, the mean equinox is at an angle, Sip — tan" x (j/^-/xy). Se = e — e 
is the nutation in obliquity, and e is the mean obliquity (the dihedral angle between the plane of the 
ecliptic and the mean plane of the equator). “Mean” as used in this section implies that the short- 
period (T < 18.6 years) effects of nutation have been removed. Actually, the separation between 
nutation and precession is rather arbitrary, but historical. The composite rotation is: 

N = A T {e)C T {SiP)A{e) (2.44) 

( cos Sip 
— cose sin Srp 
— sine&inSip 


cose sin Sip sine sin 5t/> 

cos ecos e cos Sip + sine sine cos e sin e cos S tp — sine cose 
sine cose cos Sip — cose sine sin e sin e cos S ip 4- cose cose 
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The 1980 IAU nutation model (Seidelmann, 1982, and Kaplan, 1981) is used for obtaining the 
values for Stp and e — e. The mean obliquity is obtained from Lieske et al. (1977) or from Kaplan 
(1981): 

e = 23° 26' 21."448 - 46."8150 T - 5."9 X 10 -4 T 2 + l."813 X 10~ 3 T 3 (2.45) 

This nutation in longitude (5^>) and in obliquity ( Se = e — e ) can be represented by a series expansion 
of the sines and cosines of linear combinations of five fundamental arguments. These are (Kaplan, 
1981, Cannon, 1981): 

1. The mean anomaly of the Moon: 


ai =l= 485866”. 733 + (1325 r + 715922".633) T 

+ 31". 310 T 2 + 0".064T 3 (2.46) 

2. The mean anomaly of the Sun: 

ct 2 = l' = 1287099". 804 + (99 r + 129258l".224) T 

- 0".577 T 2 - 0."012 T 3 (2.47) 

3. The mean argument of latitude of the Moon: 

a 3 = F = 335778". 877 + (1342 r + 295263".137) T 

- 13". 257 T 2 + 0".011 T 3 (2.48) 

4. The mean elongation of the Moon from the Sun: 


a 4 = D = 107226l".307 + (l236 r + 110560l".328) T 

- 6". 891 T 2 + 0".019 T 3 (2.49) 

5. The mean longitude of the ascending lunar node: 


as = O = 450160". 280 - (5 r + 482890".539) T 

+ 7".455 T 2 + 0".008 T 3 (2.50) 

where l r = 360° = 1296000". 

With these fundamental arguments, the nutation quantities then can be represented by 


and 


N 

6 '( , = J2 

y=ii 

N 

Se = ^2 

y=i 


(Aoj + A\jT) sin 


un | 

(Boj + BijT) cos (T) j j 


(2.51) 

(2.52) 


where the various values of a,-, k }i , A } , and Bj are listed in Table II. 

An additional set of terms can be optionally added to the nutations 6t/i and Se in Eqs. (2.51) and 

(2.52). These include the out-of-phase nutations and the free-core nutations (Yoder, 1983) with period 
03 f (nominally 460 days). The “nutation tweaks* Arf> and Ac are arbitrary constant increments of 
the nutation angles 6i/> and Se. Unlike the usual nutation expressions, they have no time dependence. 
The out-of-phase nutations, which are not included in the IAU 1980 nutation series, are identical to 
Eqs. (2.51) and (2.52), with the replacements sin «-► cos: 


AT 

r 6 ll 

^ = E 

(A 2l + A 3j T) cos ^ k } iOti(T) 

y=i 

S=1 -U 

N 

r 6 il 

Se = J2 

(B 2 j + B 3 jT) sin E 

II 

s=i JJ 


(2.53) 


(2.54) 
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Table II 


1980 IAU Theory of Nutation 


Index 

j 

Period 

( days ) 


B 

k } 3 


B 

Aoj A\j 

( 0 ". 0001 ) 

Bq, B\j 

( 0 ". 0001 ) 

1 

6798.4 

0 

0 

0 

0 

1 

-171996 

- 174.2 

92025 

8.9 

2 

3399.2 

0 

0 

0 

0 

2 

2062 

0.2 

-895 

0.5 

3 

1305.5 

-2 

0 

2 

0 

1 

46 

0.0 

-24 

0.0 

4 

1095.2 

2 

0 

-2 

0 

0 

11 

0.0 

0 

0.0 

5 

1615.7 

-2 

0 

2 

0 

2 

-3 

0.0 

1 

0.0 

6 

3232.9 

1 

-1 

0 

-1 

0 

-3 

0.0 

0 

0.0 

7 

6786.3 

0 

-2 

2 

-2 

1 

-2 

0.0 

1 

0.0 

8 

943.2 

2 

0 

-2 

0 

1 

1 

0.0 

0 

0.0 

9 

182.6 

0 

0 

2 

-2 

2 

-13187 

- 1.6 

5736 

- 3.1 

10 

365.3 

0 

1 

0 

0 

0 

1426 

- 3.4 

54 

- 0.1 

11 

121.7 

0 

1 

2 

-2 

2 

-517 

1.2 

224 

- 0.6 

12 

365.2 

0 

-1 

2 

-2 

2 

217 

- 0.5 

-95 

0.3 

13 

177.8 

0 

0 

2 

-2 

1 

129 

0.1 

i 

-j 

o 

0.0 

14 

205.9 

2 

0 

0 

-2 

0 

48 

0.0 

1 

0.0 

15 

173.3 

0 

0 

2 

-2 

0 

-22 

0.0 

0 

0.0 

16 

182.6 

0 

2 

0 

0 

0 

17 

- 0.1 

0 

0.0 

17 

386.0 

0 

1 

0 

0 

1 

-15 

0.0 

9 

0.0 

18 

91.3 

0 

2 

2 

-2 

2 

-16 

0.1 

7 

0.0 

19 

346.6 

0 

-1 

0 

0 

1 

-12 

0.0 

6 

0.0 

20 

199.8 

-2 

0 

0 

2 

1 

-6 

0.0 

3 

0.0 

21 

346.6 

0 

-1 

2 

-2 

1 

-5 

0.0 

3 

0.0 

22 

212.3 

2 

0 

0 

-2 

1 

4 

0.0 

-2 

0.0 

23 

119.6 

0 

1 

2 

-2 

1 

4 

0.0 

-2 

0.0 

24 

411.8 

I 

0 

0 

-1 

0 

-4 

0.0 

0 

0.0 

25 

131.7 

2 

1 

0 

-2 

0 

1 

0.0 

0 

0.0 

26 

169.0 

0 

0 

-2 

2 

1 

1 

0.0 

0 

0.0 

27 

329.8 

0 

1 

-2 

2 

0 

-1 

0.0 

0 

0.0 

28 

409.2 

0 

1 

0 

0 

2 

1 

0.0 

0 

0.0 

29 

388.3 

-1 

0 

0 

1 

1 

1 

0.0 

0 

0.0 

30 

117.5 

0 

1 

2 

-2 

0 

-1 

0.0 

0 

0.0 

31 

13.7 

0 

0 

2 

0 

2 

-2274 

- 0.2 

977 

- 0.5 

32 

27.6 

1 

0 

0 

0 

0 

712 

0.1 

-7 

0.0 

33 

13.6 

0 

0 

2 

0 

1 

-386 

1 

o 

200 

0.0 

34 

9.1 

1 

0 

2 

0 

2 

-301 

0.0 

129 

- 0.1 

35 

31.8 

1 

0 

0 

-2 

0 

-158 

0.0 

-1 

0.0 

36 

27.1 

-1 

0 

2 

0 

2 

123 

0.0 

-53 

0.0 

37 

14.8 

0 

0 

0 

2 

0 

63 

0.0 

-2 

0.0 

38 

27.7 

1 

0 

0 

0 

1 

63 

0.1 

-33 

0.0 

39 

27.4 

-1 

0 

0 

0 

1 

-58 

- 0.1 

32 

0.0 

40 

9.6 

-1 

0 

2 

2 

2 

-59 

0.0 

26 

0.0 
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Table II cont. 


1980 IAU Theory of Nutation 
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Table II cont. 


1980 IAU Theory of Nutation 


Index 

j 

Period 

(days) 

k i 1 

kj2 

kj3 

kj4 

^j’6 

(0".000l) 

Boj Bij 

(0".0001) 

81 

14.6 

-2 

0 

2 

2 

2 

1 

0.0 

-1 

0.0 

82 

5.8 

-1 

0 

2 

4 

2 

-2 

0.0 

1 

0.0 

83 

15.9 

2 

0 

0 

-4 

0 

-1 

0.0 

0 

0.0 

84 

22.5 

1 

1 

2 

-2 

2 

1 

0.0 

-1 

0.0 

85 

5.6 

1 

0 

2 

2 

1 

-1 

0.0 

1 

0.0 

86 

7.3 

-2 

0 

2 

4 

2 

-1 

0.0 

1 

0.0 

87 

9.1 

-1 

0 

4 

0 

2 

1 

0.0 

0 

0.0 

88 

29.3 

1 

-1 

0 

-2 

0 

1 

0.0 

0 

0.0 

89 

12.8 

2 

0 

2 

-2 

1 

1 

0.0 

-1 

0.0 

90 

4.7 

2 

0 

2 

2 

2 

-1 

0.0 

0 

0.0 

91 

9.6 

1 

0 

0 

2 

1 

-1 

0.0 

0 

0.0 

92 

12.7 

0 

0 

4 

-2 

2 

1 

0.0 

0 

0.0 

93 

8.7 

3 

0 

2 

-2 

2 

1 

0.0 

0 

0.0 

94 

23.8 

1 

0 

2 

-2 

0 

-1 

0.0 

0 

0.0 

95 

13.1 

0 

1 

2 

0 

1 

1 

0.0 

0 

0.0 

96 

35.0 

-1 

-1 

0 

2 

1 

1 

0.0 

0 

0.0 

97 

13.6 

0 

0 

-2 

0 

1 

-1 

0.0 

0 

0.0 

98 

25.4 

0 

0 

2 

-1 

2 

-1 

0.0 

0 

0.0 

99 

14.2 

0 

1 

0 

2 

0 

-1 

0.0 

0 

0.0 

100 

9.5 

1 

0 

-2 

-2 

0 

-1 

0.0 

0 

0.0 

101 

14.2 

0 

-1 

2 

0 

1 

-1 

0.0 

0 

0.0 

102 

34.7 

1 

1 

0 

-2 

1 

-1 

0.0 

0 

0.0 

103 

32.8 

1 

0 

-2 

2 

0 

-1 

0.0 

0 

0.0 

104 

7.1 

2 

0 

0 

2 

0 

1 

0.0 

0 

0.0 

105 

4.8 

0 

0 

2 

4 

2 

-1 

0.0 

0 

0.0 

106 

27.3 

0 

1 

0 

1 

0 

1 

0.0 

0 

0.0 


Expressions similar to these are adopted for the free-core nutations: 

Sxp = (Aqq + A 10 T) sinfu^T) + (A 20 A 30 T) cosfw/T) (2.55) 

and 

6e = (j9 0 o + B10T) cos(o>/T) + (B20 + Bs 0 T) sin(a ;/T) (2.56) 

The nutation model thus contains a total of 856 parameters: (i=0,3; ; = 1,106) and Bij (t=0,3; 

1,106) plus the free-nutation amplitudes A i0 (i=0,3), B t0 (t=0,3). The only nonzero a priori 
amplitudes are the A 0j) A ij} B oj) B iy (/= 1,106) of the 1980 IAU nutation series. 

The nutation tweaks are just constant additive factors to the angles Sip and 6e: 


Sip — * Sip + A ip 


and 


Se — ► Se 4- Ae 


(2.57) 

(2.58) 
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It is emphasized that, for the present, the default nutation model in GPSOMC is just the 1980 IAU 
nutation model. 

2.6 PRECESSION 

The next transformation in going from the Earth-fixed frame to the geocentric inertial frame is 
the rotation P . This is the precession transformation from mean equatorial coordinates of date to the 
equatorial coordinates of the reference epoch (c.g., J 2000 ). It is the transpose of the matrix given on 


page 7 of Melbourne et al. (1968): 

Pxx — cos cos ©a cos Za — sin $a sin Za (2.59) 

Pi 2 = cos (a cos ©a sin Za + sin $a cos Za (2.60) 

P 13 — cos $A sin ©A (2.61) 

P21 = — sin $a cos ©a cos Za — cos £a sin Za (2.62) 

P22 — — sin ^a cos ©a sin Za + cos £a cos Za (2.63) 

P23 — — sin ^a sin ©a (2.64) 

P31 = — sin ©a cos Za (2.65) 

P32 = — sin ©a sin Za ( 2 . 66 ) 

P33 = cos ©a (2.67) 

With the angular units in arc seconds, the arguments are 

£A = 2306.2181 T + 0.30188 T 2 + 0.017998 T 3 ( 2 . 68 ) 

Z A = 2306.2181 T + 1.09468 T 2 + 0.018203 T 3 (2.69) 

©a = 2004.3109 T - 0.42665 T 2 - 0.041833 T 3 (2.70) 


These expressions are given by Lieske et al. (1977) and by Kaplan (1981). This completes the standard 
model for the orientation of the Earth. 


2.7 PERTURBATION ROTATION 

The standard model for the rotation of the Earth as a whole may need a small incremental 
rotation about any one of the resulting axes, for example, in compensating for defects in the a priori 
precession model. Define this perturbation rotation matrix as 

O = AjAyA, (2.71) 

where 

(10 0 \ 

A a = 0 1 Se x \ (2.72) 

Vo -se x 1 J 

with SQ X being a small angle rotation about the x axis, in the sense of carrying y into s; 

f 1 0 -se y \ 

A„ = 0 1 0 (2.73) 

V$e„ 0 1 J 

with SQy being a small angle rotation about the y axis, in the sense of carrying z into x; and 

( 1 se M 0\ 

A, = -50, 1 0 (2.74) 

Vo 0 1 ) 
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with 60, being a small angle rotation about the > axis, in the sense of carrying x into y. For angles 
of the order of 1 arc second we can neglect terms of order SO 2 Re as they give effects on the order 
0.015 cm. Thus, in that approximation 

60, -60 v \ 

1 60, (2.75) 

-se, 1 ) 

60< o + 60<T + fi[T) (2.76) 

which is the sum of an offset, a time-linear rate, and some higher-order or oscillatory terms. Currently, 
only the offset and linear rate are implemented. In particular, 60 y is equivalent to a change in the 
precession constant. Setting 

60, = 60„ = 60, = 0 (2.77) 

gives the effect of applying only the standard rotation matrices. 

Starting with the Earth-fixed vector, r Eoi we have shown in Sections 2. 1-2.6 above how we obtain 
the same vector, r/ } expressed in the geocentric inertial frame: 

r/ = QPNUXY(te 0 H- A) (2.78) 


In general, 


n = -se M 
\ se v 


6 Si = SBiit) = 


2.8 GEOCENTER OFFSET AND COORDINATE SCALING 

The Earth-fixed reference frame is essentially derived from VLBI measurements, which are com- 
pletely insensitive to the location of the Earth’s center of mass. Practically, receiver coordinates are 
based on the location of a reference station, which is derived either from spacecraft tracking data or 
satellite laser ranging (SLR) measurements. Consequently, there might be sisable errors in fiducial 
station locations (~ 10 m or ~ 1 m for the two techniques, respectively). GPS tracking data is capable 
of determining the position of the geocenter to an accuracy dependent on the quality and quantity of 
the range observations. In order to allow for the possibility of solving for the location of the geocenter, 
Cartesian offsets were introduced, such that the receiver coordinates (including those of the reference 
station (s)) are 


x Eo X E 0 + X GG 

VEo — * yEo + VGC (2.79) 

ZEo * z Eo + Zqc 

or in vector notation, 

te 0 r E 0 + *gc (2.80) 

Another problem with coordinate systems at the present level of understanding is the uncer- 
tainty that the coordinate transformations for a priori satellite orbits are being carried out correctly. 
Introduction of a scale factor a, which multiplies Earth-fixed coordinates relative to the geocentric 
inertial coordinates of the satellites, permits empirical detection of such problems. It gives rise to the 
transform at ions 


or 


X E 0 “ * ax E 0 


yE 0 a VEo 

(2.81) 

ZEo az E 0 


TEo —* aT E 0 

(2.82) 
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in vector notation, for the coordinates of the receivers. 

Upon inclusion of both of these modifications in the Earth-fixed station locations r^ 0 of Eq. (2.78), 
the final expression for the station location vector in geocentric inertial coordinates becomes 

17 = nPNUXY(otT Eo + tgc + A) = Q(ar Eo + t G c + A) (2.83) 


2.9 PHASE CENTER OFFSETS 

Geometric offsets may exist between transmitter and receiver phase centers and the standard 
reference points. In the case of the ground-based receiver, this offset takes the form of a “site vector* 
from a surveyor’s benchmark to the receiver phase center. Normally, the position of the benchmark, 
Ybm j is desired, rather than the position of the phase center, rp G . These two vectors are related 
through the site vector tsurv* 


*E 0 = = Tpo + TSUJSV (2.84) 

For GPS satellites, the ephemeris reference point is the spacecraft center of mass (CM). The vector 
offset between the phase center and CM is given by 0.211 i + 0.886 k meters in a coordinate system in 
which the unit vector k points from the spacecraft CM to the Earth center, j is the normalized cross 
product of k with the unit vector from the spacecraft to the Sun, and i completes a right-handed 
system (Winn, 1984). Both of these geometric offsets are incorporated into GPSOMC modeling. 
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SECTION 3 


OBSERVABLES AND CLOCK PARAMETERS 


Range measurements from GPS satellites are based on the detection of the phase of transmitted 
electromagnetic signals. The physical observable is the difference between the reference phase of the 
transmitter at the time of signal emission and the reference phase of the receiver at the time of signal 
reception. Such measurements of necessity involve detailed and careful consideration of a number of 
time scales. These are intimately related to the definitions of the observables, and for that reason this 
section contains detailed descriptions of both the observables and clock models. 

The GPS satellite transmissions are one-way transmissions originating at the spacecraft. A carrier 
signal modulated by a pseudo-random noise code is broadcast (Spilker, 1978). Both a “carrier phase” 
observable <p and a “pseudo-range” observable R are modeled in GPSOMC. The observed observables are 
related in a prescribed way to the physically detected phase. The computed observables are defined 
in terms of clock differences since the reference phase at the receiver is derived from the station clock 
and the reference phase at the transmitter is derived from the space vehicle clock. For both carrier 
phase and pseudo-range the computed observable is given by the theoretical difference between the 
space vehicle and station clocks plus a bias term. The clock difference is accumulated as the sum of 
five components — the first is the difference between station clock time and proper time at the station; 
the second is the difference between proper time and coordinate time at the station; the third is the 
difference between the coordinate time of signal emission and the coordinate time of signal reception; 
the fourth is the difference between coordinate time and proper time at the spacecraft; the fifth is the 
difference between proper time at the spacecraft and spacecraft clock time. 

These terms are readily interpreted. The first term is referred to as a station clock error and 
the fifth term is referred to as a space vehicle clock error. The second and fourth terms are general 
relativistic time transformations. The third term is the solution to the light time equation connecting 
the events of signal emission and signal reception. Of these five contributions the third term is 
normally dominant. This term is also referred to as the geometric range between spacecraft and 
receiver, whence the use of “range” in referring to GPS observables. Tropospheric and ionospheric 
delays are not considered here; they are discussed in Sections 1 and 4 of this report. The observables 
are assumed to have been calibrated to remove the effects of instrumental delays. A detailed model 
for instrumental phase delay is not provided, but a bias term is included in the definition of both data 
types. This bias might account for 

(i) an unknown integer number of cycles of phase, 

(ii) an uncalibrated offset in the absolute phase of an oscillator, 

(iii) uncalibrated delays within transmitter or receiver electronics, or 

(iv) an uncalibrated offset between the phase reference used for signal detection and the phase 
reference used for time-tag generation within a receiver. 

The geometric portion of the observable is calculated in the geocentric reference frame, moving 
with the Earth, with inertial J2000 oriented axes. This reference frame is the generalization of the 
local inertial reference frame containing the Earth (Ashby and Bertotti, 1984) and provides a simple 
setting for applying all necessary relativistic corrections. Terrestrial Dynamic Time (TDT) is used as 
coordinate time (Thomas, 1975 and Ashby and Allan, 1979). The spacecraft equations of motion are 
integrated in this reference frame using TDT as the independent variable. As described in Section 2 
of this report, station coordinates are transformed into this reference system. Care must be taken 
when selecting a priori values for model parameters since particle mass and physical distance in 
the geocentric system differ by a scale factor from current determinations of same in the celestial 
coordinate system (Hellings, 1986). 

We use *2 to denote the time of signal emission and *3 to denote the time of signal reception. 
Let £3, £3, and (3 refer to, respectively, station clock time, proper time, and coordinate time at the 
receiver. Let *2, *2, and *2 refer to, respectively, space vehicle clock time, proper time, and coordinate 
time at the transmitter. 
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5.1 PSEUDO-RANGE 


Observed pseudo-range is given by 


Z° = c ( f 3 - t 2 ) (3.1) 

with time tag t 3 , where the times and are the actual times kept by the station and spacecraft 
clocks, respectively. Correlation of the received pseudo-random noise code generated and transmitted 
by the space vehicle with a local copy of the code within the receiver allows for direct conversion 
between the detected received phase and clock readings. 

Computed pseudo-range is given by 

Z G = c ( t 3 - t 3 + Bstn ) (3.2) 

with time tag t 3 , where now the times t 3 and t 3 are the modeled times of, respectively, the station 
and space vehicle clocks. To facilitate computation, this expression is rewritten as a sum of six terms, 

Z c — c(t 3 -h 
+ t 3 — t 3 
4- t 3 — t 3 
4- t 3 — t 3 
4 - t 3 — 1 3 

4- Bstn ) (3.3) 

The first five terms are the five clock differences discussed above, and the last term is the bias discussed 
above. 

The station clock error is modeled as a quadratic function, 

t 3 -h = asTN'Z + &STJV,*(*3 *3 0 ) + C STN,Z (*3 “ *3 0 ) (3.4) 

On the right-hand side t 3o is a specified epoch. The coefficients <istn,r> r, and cstn,r are 

generally unknown but presumed to be small. Ideally, in geocentric coordinates, all clocks fixed to 
the surface of the Earth run at approximately the same rate, with the precise rate depending weakly 
on the geodetic coordinates of the clock. We assume that all clocks are adjusted as necessary so that 
their average rate is the same as that of coordinate time. The rate of coordinate time has been defined 
to agree with the SI second (Kaplan, 1981). Thus, in our model, station proper time, ideal UTC time, 
and coordinate time all run at the same rate. This is not a limitation since the deviation of the rate 
of any clock from the assumed rate can be modeled as a station clock error. 

The difference between proper time and coordinate time at the station is given by 

f 3 - h = - [ {TAI - UTC) + {TDT - TAI) ] (3.5) 

where TAI — UTC is an integer number of leap seconds which changes approximately once a year 
and TDT — TAI is defined to be 32.184 sec. 

The geometric range is given by the solution to the light time equation, 

. „ _ I rsTu(h) -rsv(ta)| , Ax 

* 3 ~ *3 = + A *r el ltom (3.6) 

where rsrtf(t 3 ) i« the inertial receiver position at the time of signal reception, rsy(£ 2 ) is the space 
vehicle position at the time of signal transmission, and A t re \ gtcm is the general relativistic correction to 
the geometric range. To begin the computation of the range, we start from the measurement epoch f 3 . 
This is converted to t 3 using Eq. (3.4) and to t 3 using Eq. (3.5). The station position in geocentric 
inertial coordinates is then calculated at the epoch t 3 using Eq. (2.83). Thus station clock errors 


18 



affect the computation of the geometric range. Given the station position and a PV file containing 
the geocentric inertial spacecraft position as a function of coordinate time, the light time equation is 
iteratively solved to obtain 1 2 - 


/*+*) 

c 2 




*3 -4 f) - r !a/ e_ At ri. 


1 — r 


(0 

12 


*V(^) 


(3.7) 


The space vehicle position used in this computation is that of the antenna phase center, as discussed in 
section 2.9 of this report. The gravitational delay A t re i g , om is given by (Tausner, 1966 and Holdridge, 


1967) 


At re l s 


«om 


(1 + 1ppn)pe ^ ri + + ri 2 

c 3 r x + r 2 - r i2 


(3.8) 


where Ippn is 1 for general relativity, ^ is the Earth’s gravitational constant, and, as in Eq. 


(3.7), 


*1 = l*5rjv(t3) 
*2 = rsv(*2) 

r i2 = *i -t 2 


n = ki| 
*2 = 1 
*"i2 — |ria| 


(3.9) 


We use r 5 v(< 3 ) as the initial estimate for rsv(* 2 ); the light time iteration converges very fast. 

The rate of an ideal space vehicle clock will differ from the rate of coordinate time due to the 
difference in gravitational potential and speed between the space vehicle clock and the reference clock. 
This rate difference may be decomposed into two components: a bias component which depends only 
on the nominal semi-major orbit axis of the space vehicle, and a periodic component which depends 
mainly on the orbit eccentricity. The frequencies of the GPS satellite clocks were purposely set 
low before launch, relative to the nominal published values, to offset the nominal bias rate difference 
between space vehicle time and coordinate time (Spilker, 1978). Accordingly, our model rate difference 
between proper time and coordinate time at the spacecraft contains only a periodic component. This 
component depends only on orbit eccentricity. A bias component due to off-nominal semi-major orbit 
axis and a small periodic component due to Earth oblateness and solar perturbations are absorbed in 
the space vehicle clock error model. The difference between proper time and coordinate time at the 
space vehicle is given by 


t 2 - h = (TAI - UTC) + (TDT - TAI) + A t relclk (3.10) 

where the clock correction term is 


At r el clk = 2 r5v(*2) (* 2 ) / C* 


(3.11) 


This result may be derived from integration of the usual differential equation relating coordinate 
and proper time [e.p., Eq. (l) of Thomas (1975)] in the geocentric reference system, modeling the 
space vehicle orbit about a point-mass Earth as elliptical while ignoring the differential gravitational 
potential due to perturbing bodies. The resulting rate expression is accurate to about one part in 10 12 
for the GPS orbits, which is comparable to the instability of the GPS satellite clocks. Note that the 
large (» 50 sec) defined offset between coordinate and proper time appears in Eqs. (3.5) and (3.10) 
with opposite sign, and hence does not affect the computation. It is necessary, though, to include this 
offset due to our convention of basing observable time tags on UTC while indexing the PV file by 
TDT. 

Analogously to the station clock error, the space vehicle clock error is also modeled as a quadratic 
function, 

h - h = + t>sv t z{i2 — *2 0 ) + —t 2 0 ) 2 ( 3 * 12 ) 
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On the right-hand side t 2o is a specified epoch. The space vehicle clock error does not affect geometric 
range. The last term in the model is the pseudo-range bias Bstn » which depends on the station but 
is independent of the spacecraft. It also does not affect the value of the computed geometric range. 

3.2 CARRIER PHASE 

Observed carrier phase is given by 

<p° = (— c/w n ) (<f>sv — 4 >stn) (3.13) 

with time tag t 3 , where <f> S v is the received RF (radio frequency) phase of the space vehicle signal 
at time £ 3 , <f>sTN i® the receiver reference phase at time t 3 , and w n is the nominal value for both the 
transmitted frequency of the space vehicle signal and the receiver mixing frequency. 

The receiver reference phase is modeled as 

<t>STN{t3) = ( f 3 - t' 3o ) (3.14) 

Since phase is a physical invariant, the received phase of the space vehicle signal at the time of signal 
reception is the same as the space vehicle reference phase at the time of signal transmission, 

<fisv(h) = <f>ref(t2) (3.15) 

The space vehicle reference phase is modeled as 

^rc/(t 2 ) — ( ?2 — ?2o) (3.16) 

Thus the computed carrier phase is given by 

<p° = c ( ? 3 - £3 + Rsv^rjv) (3.17) 

i 

with time tag f 3 . The integration constants and have been absorbed in the carrier phase bias 
Bsv,stn • Deviation of the station mixing frequency from nominal is modeled as a station clock error, 
and deviation of the space vehicle transmitter frequency from nominal is modeled as a space vehicle 
clock error. Note that the modeled value for carrier phase appears to be the same as the modeled 
value for pseudo-range, except for the substitution of the carrier phase bias for the pseudo-range bias. 
A further distinction is that we allow the modeled values of station time and spacecraft time to be 
different for the two data types. 

Just as for pseudo-range the expression for computed carrier phase is rewritten as a sum of six 
terms, 


<P C = C ( t 3 - t 3 


+ £3 — t 3 

+ “ t 2 

H~ t 2 — £2 
H- ~ £ 2 

+ Bsv,STN ) 


( 3 . 18 ) 


The station clock error is modeled as 


t 3 -t 3 = a S TN,ifi + f>STN,v(t3-t 3 o) + C S TN,v (<3 ~ t 3o ) 2 


( 3 . 19 ) 
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and the space vehicle clock error is modeled as 

t 2 - t 2 — cLsv.ip + bsv t p{t 2 — ho) + csVtpih - ho) 2 (3.20) 

The carrier phase bias Bsv.stn depends on both the station and the space vehicle. The other terms 
in Eq. (3.18) are evaluated just as for pseudo-range. They will, however, have different numerical 
values unless the station clock error coefficients are identical for the two data types. The value of the 
station clock error affects the geometric range, while the values of the carrier phase bias and space 
vehicle clock error do not. 

In cases where the behavior of transmitter and receiver clocks permits, both the range and phase 
observables may be modeled with a common set of parameters. By analogy with Eqs. (3.4), (3.12), 
(3.19), and (3.20), the common clock errors may be written as 

£3 — = a STN + &STAr(£3 — tSo) + C STN (*3 ~ *3 0 ) 2 (3.2 1) 

and 

*2 ~ h = asv + &sv(t2 — f2o) + c sv (fa “ ^2 0 ) 2 (3.22) 
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SECTION 4 


TRO P O SPHERE 

The delay experienced by the incoming signal due to the Earth’s atmosphere can be modeled 
using a spherical-shell troposphere consisting of a wet component and a dry component. If E is the 
apparent geodetic elevation angle of the spacecraft, the tropospheric contribution to the range is 

Ptrap - PZt rv Rdry{E) + PZ W9t ^wet {E) (4.1) 

where pz is the additional delay at local senith due to the presence of the troposphere, and R is an 
elevation angle mapping function. For recent VLBI data, it was found that modeling the zenith delay 
as a linear function of time improves troposphere modeling considerably. The dry and wet zenith 
parameters are therefore written as 

PZi.w = P%t.w + p**.«(* " *o) ( 4 - 2 ) 


where to i* a reference time. 

The analytic mapping function developed by Lanyi (1984) is used for mapping zenith values to 
the observed elevation angles. In its simplest form, this mapping function employs average values of 
atmospheric constants. Provision is made for specifying surface meteorological data acquired at the 
time of the experiments, which may override the average values. An approximate partial derivative is 
obtained with respect to one parameter in the Lanyi mapping function; this permits adjustment even 
in the absence of surface data. 

Here we attempt to give a minimal summary of the final formulas. The tropospheric delay is 
written as: 

Puop = F(E)/siaE (4.3) 

where 


F(E) = PZ iry Fdry(E) + pZ w .,Fwet{E) 

+ [/>£<„ + 2pz irt Pz„ME) + 4.* r„(£)]/A + p 3 Zirv F M (E)/ A 2 (4.4) 

The quantities pz 4rv an d Pz w0t have the usual meaning: zenith dry and wet tropospheric delays. A 
is the atmospheric scale height, A = kTo/mg ct with A: = Boltzmann’s constant, T 0 = average surface 
temperature, m = mean molecular mass of dry air, and g c = gravitational acceleration at the center of 
gravity of the air column. With the standard values k = 1.38066 x 10" 16 erg/K, m = 4.8097 x 10" 23 
g, g c = 978.37 erg/gcm, and the average temperature for DSN stations To = 292 K, the scale height 
A = 8567 m. 

The dry, wet, and bending contributions to the delay, F<i rv (E) t F wet (E) y and Fm, 62 . 63,64 {E) y are 
expressed in terms of moments of the refractivity as 


Edry(E) — Aio^jG^AAfno, u) + 3auM 2 ioG r3 (Afiio, u)/4 (4.5) 

F wct (E) = Aoi(E)G(\Mioi/Moouu)/Moox (4.6) 

Fbi{E) = [crG 3 (Muo, txj/sin 2 £ — Afo2o^ 3 (Afi2o/^020j tx)]/2 tan 2 E (4.7) 

Fb2(E) = ”MonG 3 (Miii/Afoii, u)/2Afooi tan 2 E (4.8) 

•F* 63 {E} — -M 002 G 3 (M X 02 /Afoo 2 , u)/ 2 Af 0 2 01 tan 2 E (4.9) 

F 64 (-^) = Mo 3 oG 3 (A/ 13 o/M 0 3 o, u )/2 tan 4 E (4-10) 

A misprinted sign in the last of Eqs. (5) of Appendix B of Lanyi (1984) has been corrected in 

Eq. (4.10). Here G(<?,u) is a geometric factor given by 

G{q,u) = (l + gu) -1/2 (4.11) 
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with 


u = 2a/ tan 2 E (4.12) 

where a is a measure of the curvature of the Earth’s surface A jR with standard value 0.001345. 

The quantities Ai m (E) and Mum are related to moments of the atmospheric refractivity, and are 
defined below. Aio(£ r ) involves the dry refractivity, while Aoi(^) is the corresponding wet quantity. 
The Ai^E) are given by 


4 (F\ = M , Y>A (-ir +fc (2n-l)!!Af n - Mi 

A tm {E) Moim + 2_^2^ 2 n k\[r\ — A:)! 


n= 1 fc= 0 


■J 

U 

n 

XMilm 


1 + XuMum/Moim 


M 0 Im 




(4.13) 


with the scale factor A = 3 for E < 10° and A= 1 for E > 10°. Only the two combinations (/, m) = 
(0,1) and (1,0) are needed for the Ai m (E). The moments of the dry and wet ref r activities are defined 
as 

CO 

Vniy = J q n f' dry {q)fl )e t(q) d <l (414) 

0 

where fdry t u/et(<7) are the surface-normalized refr activities. Here, n ranges from 0 to 1, i from 0 to 3, 
and j from 0 to 2; not all combinations are needed. Carrying out the integration in Eq. (4.14) for a 
three-section temperature profile gives an expression for the general moment Af n jy: 

M„ <y /n! = (1 - e~ cu,l )/a n+1 + [l - 7? + " +1 (?i,<7 2 )] f[ j-f-y 

» = 0 

+ ( -^T* +n+1 ( qi ,q 2 )/a n+ 1 (4.15) 


Here, 


$ 2 ) = 1 — (<?2 — ?x)/ Qt 


(4.16) 


The quantities qi and q 2 are the scale-height normalised inversion and tropopause altitudes, respec- 
tively. For the standard atmospheric model, qi = 0.1459 and q 2 = 1.424. The constants a and b are 
functions of the dry (a = 5.0) and wet (fi — 3.5) model parameters, as well as of the powers of the 
refractivities (t and j) in the moment definitions. Table III gives the necessary a’s and 6’s. 


Table III 

Dependence of the Constants a and 6 
on Tropospheric Model Parameters 


t 

j 

a 

b 

1 

0 

1 

a — 1 

0 

1 

0 

afi- 2 

2 

0 

2 

2(a-l) 

1 

1 

p+i 

/9(a + 1) - 3 

0 

2 

2/9 

2(a/9 - 2) 

3 

0 

3 

3(a-l) 


Note that the normalization is such that Af 0 xo = 1; this moment has therefore not been explicitly 
written in Eqs. (4.5)— (4.10). 
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At present, provision is made for input of four meteorological parameters to override the standard 
(average) values of the Lanyi model. These are (1) the surface temperature To, which determines the 
atmosphere scale height (default value 292 K); (2) the temperature lapse rate W , which determines 
the dry model parameter a (default values W = 6.8165 K/km, a = 5.0); (3) the inversion altitude 
hi, which determines q± = h\j A (default value hi — 1.25 km); and (4) the tropopause altitude h 2 , 
which determines = h 2 /A (default value h 2 — 12.2 km). A fifth parameter, the surface pressure 
poi is not used at present. Approximate sensitivity of the tropospheric delay (at 15° elevation) to the 
meteorological parameters is —0.3 cm/K for surface temperature, 1 cm/(K/km) for lapse rate, and 
— 1 cm/km for inversion and 0.2 cm/km for tropopause altitude, respectively. 
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SECTION 5 


DERIVATIVES OF OBSERVABLES WITH RESPECT TO MODEL PARAMETERS 


This section presents expressions for the partial derivatives of the computed observables with 
respect to the three classes of model parameters described in this document, as well as with respect 
to dynamic parameters related to the spacecraft state. The relativistic corrections, the tropospheric 
correction, and the space vehicle clock error term in the computed observable depend weakly on 
geometric parameters, on the station clock error, and on dynamic parameters. These functional 
dependences are neglected when computing partial derivatives. 

We begin with some notation. From Eqs. (3.2) through (3.17) the computed observable can be 
written as 

C == c ( £3 — £ 3 ) + p + c ( 1 2 — fa) + c B v (5*1) 

where C is Z c for pseudo-range and <p c for carrier phase, and 1 / is STN for pseudo-range and 
SV y STN for carrier phase. 

The clock error terms are given by 

^3 ^3 = a STN„p + &STJV,0(*3 ~ *3 0 ) + CSTN,0 (*3 — *3 0 ) 

and 

?2 ” h r = <*sv,$ + bsv.pih — ?2 0 ) + c svrf{t 2 — ^2o) 3 (5.3) 

where fi is R for pseudo-range and <p for carrier phase. The geometric range is given by 

P = I r S TN - Tsv I (5.4) 

where Tstn is the geocentric inertial station position at the time £3 of signal reception and rsv is 
the space vehicle position, in the same system, at the time £3 of signal emission. The station position 
comes from Eq. (2.83), 

r STN — CIPNUXY ( OiTsTNo + TqC + A) (5*5) 

The space vehicle position and the times of signal emission and reception are related through the 
light time equation, from which we may write the geometric range as 

p = c ( £3 — £3 ) 

and, equivalently, 

c 2 (*3 - *2) 2 = ( tstn — rsv) * ( tstn — *sv) (5.7) 



(5.2) 


5.1 GEOMETRIC PARAMETERS 


The class of “geometric* parameters includes receiver coordinates and their rates, Love numbers 
and tide phase, relativistic 7 , geocenter offsets, coordinate scale factor, UTPM parameters, nutation 
amplitudes, and the coefficients of the perturbation rotation. Computed observables for both the 
carrier phase and pseudo-range have identical functional dependence on these geometric parameters 
through the geometric range. The partials will therefore be written for the geometric range given by 
Eq. (5.4). Symbolizing one of the geometric parameters by ij, we have 


dp_ 

dr} 


— -( tstn — *sv) T 
P 

= -(rsTAr - rsv) T 
P 

— -{tstn - Tsv) 7 

p 


q^( T STN ~ r $v) 

[ &tstn _ 9tsv dr stn 1 

dt) dvsTN dt) J 

[j 3 _ drsv 1 dr stn 
I* dtsTN •* dr\ 


(5.8) 
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The light time correction term is computed first: 


d*sv 

&*STN 


= Tsv 


dt 2 

drsTN 


(5.9) 


The partial on the right-hand side is obtained from Eq. (5.7) by implicit differentiation: 


2 c 2 (t 3 - t 2 ) = 2 (rsrw - r S v) r [ 7 a “ *sv 


dtn 


&tstn J 


(5.10) 


Solving for 


dt 2 

3r stn 


gives 


3*2 

drsTN 


—{tstn ~ rsy ) 1 


cp 


(*- 


*STN — TSV 


TSV ' 
C - 


(5.11) 


The computation of the geometric partials is complete except for the parameter-dependent terms 
Qtstn 

dtf 


Taking the parameters in order of their appearance on the right-hand side of Eq. (5.5), the partial 
derivatives are 


Ststn 

3 

&TsTN 

3 


dn 

3 6G x , ytX 

an 

3 £0* ll/t . 


-PNU XY {ax st N 0 + tgc + A) 

o 

PNU XY (arsTNo + r ac + A) 


(5.12) 

(5.13) 


for the perturbation rotation parameters; 

^^- = n / > (|^)A-y(arsTi V ,+r G c + A) (5.14) 

?j^=np(^) XY (aT STNo +TGc + A) (5.15) 

for the nutation series amplitudes; 

^^ = np (^) xy(arsTWo+rGC + A) (516) 

= ° p ( yz) xY{arsTN ° +rac + a) (5i7) 


for the nutation offsets; 


Ststn 

d{UT 1 - UTC) 


= nM 3(Wl- UTC) XY{aTsTN ° + Toc + A) 


dH 


= OPN^-XY(ar ST N 0 + tgc + A) 
oH 


for universal time, and 


dr S TN 

301,2 


/ 3 xy \ . 

= n PNU ( dQ^ 2 ) ( ar 5 T // 0 + r ac + A) 


(5.18) 

(5.19) 


(5.20) 
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d T STN _ _ ftpffU ^ \ [ aT g TNa + XqC + A) 

ae li2 x d&i, 2 ' 


(5.21) 


for polar motion. The only quantities remaining to be evaluated in Eqs. (5.12— 5.20) are the partial 
derivatives of the various rotation matrices. These are exactly the same as the partials entering 
modeling of VLBI observations, and are set out in some detail in the publication of Sovers and 
Fanselow (1987), Sec. 2.12. For example, the first partial in Eq. (5.12) is 



(5.22) 


The partials involving UTPM rates are 


dX , 


BY , , 

de 2 ~ (t w) 


3U . , 

3H~ (t W) 




0 0 

— sin 02 cos 02 

— cos 02 — sin 02 

in H — cos H 0 \ 
a H —sin H 0 J 

0 0 0 J 


(5.23) 


(5.24) 


(5.25) 


Continuing the evaluation of partial derivatives with respect to geometric parameters, and using 
Q = QPNUXY from Eq. (2.21) to simplify notation, we have for the station coordinates and their 
time rates of change, 


For the geocenter offset components, 
For the coordinate scale factor, 


“ r 5rJV + *5 TN (* “ *o) 

(5.26) 

= Q(arsTN 0 + t gc + A) 

(5.27) 

3r STN * 

dl %TH ~ a 

(5.28) 

IP 2 - -(<-•. w 

ot STN 

(5.29) 

dr S TN q 

dr ac V 

(5.30) 

drsTN 

da = Qrs ™° 

(5.31) 


The only parameters in the tidal displacement A are the two solid- Earth- tide Love numbers h 
and /, and the tide phase angle ip. The corresponding partial derivatives of the range are somewhat 
complicated: 

Btstn 


dh 

d*STN 

dl 


— QVW 


= QVW 



(5.32) 

(5.33) 
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( dgx(Q \ 

dgaji) 

V av-i / 

where the tide phase partials of the <^s of Eqs. (2.9—2.11) are 

” R 6 ^ r P ‘ ^ ( X *yp YaXp) 


dpi. _ 3^r| 


a^» 






( r p ■ R*)( x p^» + VpF*) ~ (F«x p — X 8 y p f 


&93s _ 3/i#fp j / \r \r _ \ 

a^> - r . 6 * * x ‘ yp y * p ) 


\J x p + Vp Z, - Zp (2r p R, - * p Z.) 
\Ap + # 


(5.34) 


(5.35) 

(5.36) 

(5.37) 


Finally, the partial derivative of range with respect to the post-Newtonian parameter 7 PPN [see 
Eq. (3.8)] is 


dp _ Pe_ ln _ri +r 2 + ri 2 


#7PPAr c 2 ri + r 2 - r i2 


(5.38) 


6.2 CLOCK PARAMETERS 

Clock parameters include the coefficients of the station and space vehicle clock error models and 
also the pseudo-range and carrier phase biases. From Eq. (5.1), the partial of pseudo-range with 
respect to the first station clock error parameter is 


dZ c 


= c 


dasTN, k dasTN.k 

= c- P 


d t ~ - . dp 

(*3 — * 3 ) + 


dU 


dt$ dasTN y k 


where we have used 


dasTN, k 
dt z d 


dasTNyk dasTNyk 


(is — £ 3 ) = 1 

[t 3 — (f 3 — t 3 )] = —1 


and 


dp_ 

dt 3 


The range rate p is, from Eq. (5.6), given by 

- £) 

The partial on the right-hand side is obtained from implicit differentiation of Eq. (5.7): 
2 c 2 ( f 3 - t 2 ) (l - = 2 ( t S tn - rsv ) • (r stn - rsv 


(5.39) 

(5.40) 

(5.41) 

(5.42) 

(5.43) 

(5.44) 


28 


(5.45) 


dt 2 . 

Solving for — — gives 
ats 


dt 2 _ Cp - (rsTN - Tsv ) • TSTN 
dtz cp — (TSTN — Tsv) ' Tsv 


where we have used Gq. (5.6). Substituting into Eq. (5.43) gives 

. _ (rsrAr ~ rsv ) ■ (tstn - *sv ) 
P ~ (i _ r STN - Tsv Tsv \ 
P \ p c / 


(5.46) 


The station velocity is obtained from the time derivative of Eq. (5.5), and the space vehicle velocity 
is obtained from interpolation of the PV file. 

The partial with respect to the first space vehicle clock error parameter is 

-p— = C P~[h - h) = C (5.47) 

oa S v } R o<*sv>R 


The pseudo-range bias partial is 


dk G 

dBsTN 


(5.48) 


The pseudo-range partials with respect to the other clock error parameters are similarly computed, 
as are the carrier phase partials; here we summarise the results. 


dZ c 

d&STN,R 


dk C 

dl>STN,R 

dk c 

dCsTN,R 


(c-p)(t 3 -f 3o ) 

= (c - p)(*3 - *3 0 ) 2 


d&sv, je 

3Jg c 

dbsv,R 

aje c 

dcsv,R 


= c(t 3 - ta 0 ) 
= c(t2 — t2o ) 2 


(5.49) 


dBsTN 


(5.50) 


d<p c 

= (c - p) 

d(p c 


d&STN ,<p 

d<*sv t v 


d<f> c 

0 

1 

1 

O 

II 

d<p c 

— c(ia — ta 0 ) 

dbsr N,v 


d<p c 

II 

ft 

1 

09 

1 

O 

d<p c 

= c(ta - * 2 0 ) 

dcsTN^v 

dcsv,p 


d<p c 




dBsv,STN 


(5.51) 


(5.52) 
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5.3 TROPOSPHERE PARAMETERS 


Partials of the range with respect to the dry and wet zenith delays are obtained from Eqs. (4.3) 
and (4.4): 


dp 

dpZiry 

dp 

dpz„. t 


= [F drv {E) + 2 PZirv F bl (E)/ &}/ smE 

+ \lpz v . t Fb 2 (E)l A + Zp 2 Zirt Fbi{E) / / sin£ 

= [^toet(-E) + 2pz try Ff,2(E)/A + 2pz w „ Fbs(E) / A] / ain E 


For the senith rate parameters: 


dp 

di>Zt, m 


(t - t 0 )R d , w 


(5.53) 

(5.54) 


(5.55) 


In the analysis of data for which meteorological parameters are not available, it is convenient 
to introduce an approximation to the mapping function [Eqs. (4.3) and (4.4)] which involves a 
one-parameter estimate. This parameter p accounts for deviations from standard meteorological 
conditions. The tropospheric range is expressed as 


r\ 

P ~ (PZdry + Pz w .t )/ 8i nE + 


(5.56) 


where the partial derivative is 

dp __ [pzt rv + pz vet )uM no 

dp G(M ll0 , u) [l + G(Af 110j «)] sin E 

^ Pz* 0 t u {M no — Mxqi/Mqoi) 

G(Mno i u)G(Mioi/Mooi,u) [(7(Afno, u) + G(M 1 qi/M 0 oi 1 u)] sinl£ 


(5.57) 


5.4 SATELLITE PARAMETERS 


The final class of parameters in modeling GPS range measurements includes the dynamic pa- 
rameters characterizing the forces on the satellites, and the six “epoch state” position and velocity 
parameters, r S v 0 > *sv 0 for each space vehicle. Carrier phase and pseudo-range are affected in the same 
way, through the geometric range p, by these parameters. Here we use £ to represent any parameter 
of this class. The partial derivative of geometric range with respect to £ is calculated from Eq. (5.6): 


dp __ dt2 


(5.58) 


The partial on the right-hand side is obtained from implicit differentiation of Eq. (5.7): 

2c 2 (ta - h) = 2(rsxw - rsv) T (— ^rsv) 

= -2(r stn ~ rsv ) r | t + * sv ’^) ( 5 . 59 ) 
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dh . 

Solving for — gives 


dt 2 


(r gr ,- r ,vF^f| t> 

^ _ r STN ~ *SV *SV j 


Substituting this result into Eq. (5.58) gives 


dp 

9* 


“( *STN — 


T dr sv 


dt It, 


■(>- 


r STN — TSV 


r -f) 


(5.60) 


(5.61) 


The space vehicle velocity T 5 k and the variational partial 

t< 2 * Detailed descriptions of the spacecraft force models and 
Mathematical Description document (Wu et al. } 1986). 


are read from the PV file at time 

d£ 1*3 

variational partials are found in the OASIS 
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APPENDIX A 


Parameter Units and Physical Constants 


Early during the development of software for the GIPSY project, it was decided that the units 
of all parameters and partial derivatives shall be kks (kilometers, kilograms, seconds), instead of the 
more customary mks. With the exception of clock parameters and range biases, we adhere to the kks 
convention. To prevent numerical instabilities in parameter estimation arising from the disparity of 
magnitudes of the clock parameters and the observables, clock offsets and range biases are in units of 
jxsec (10“ 6 sec), rates in 10” 12 sec/sec, and rate rates in 10~ 18 sec/sec 2 . All lengths were measured 
in units of km, velocities in km/sec, and accelerations in km/sec 2 . Angular quantities, such as tide 
lag, UTPM, the perturbation tweaks, and nutation angles are in radians (or rad/sec for their rates). 
Finally, the senith troposphere quantities are measured in km, and their rates in km/sec. These 
conventions do not necessarily apply to the formatted GPSOMC input file, where units may be used for 
some parameters that make them more easily recognisable. 

We have tried to use the constants recommended by the IAU project MERIT (Melbourne et al., 
1983). Those that have not been defined in the text above, but which affect the results that are 
obtained using GPSOMC, are given below: 


Symbol 

Value 

Quantity 

c 

299792.458 

Velocity of light (km/sec) 

Re 

6378.140 

Equatorial radius of the Earth (km) 

oje 

7.2921151467 x 10“ 6 

Rotation rate of the Earth (rad/sec) 

f 

298.257 

Flattening factor of the geoid 

h 

0.609 

Vertical Love number 

l 

0.0852 

Horisontal Love number 

9 

980.665 

Surface acceleration due to gravity (g/cm 2 ) 
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APPENDIX B 


GPSOMC Parameters 


Table B.I 

Glossary of GPSOMC Parameters 


Parameter 

GPSOMC name 

Definition 

Reference 


SAT PEPOvvvvvv 

Coefficients in 

(3.20) 


SAT PRATvvvvvv 

space vehicle clock 

(3.20) 


SAT PRRTvvvvvv 

model (phase) 

(3.20) 

&STN,v> 

STA PEPO ss 

Coefficients in 

(3.19) 

bsTN^v 

STA PRAT 88 

station clock 

(3.19) 

CSTNyip 

STA PRRT 88 

model (phase) 

(3.19) 

Bsv,stn 

BIAS *vv_B8 ; nn 

Carrier phase bias 

(3.18) 

&SV, k 

SAT REPOvvvvvv 

Coefficients in 

(3.12) 

bsv, k 

SAT RRATvvvvvv 

space vehicle clock 

(3.12) 

Csv y z 

SAT RRRT vvvvvv 

model (range) 

(3.12) 

&STN,k 

STA REPO 88 

Coefficients in 

(3.4) 

l>STN, JZ 

STA RRAT 88 

station clock 

(3.4) 

CSTN,ft 

STA RRRT 88 

model (range) 

(3.4) 

Bstn 

BIAS PSR ss 

Pseudo-range bias 

(3.3) 

a SV,r 

SAT CEP 0 vvvvvv 

Coefficients in 

(3.22) 

bsv t <p 

SAT CRATvvwvv 

space vehicle 

(3.22) 

c SV,v 

SAT CRRTvvvvvv 

clock model 

(3.22) 

&STN,<p 

STA CEPO 88 

Coefficients in 

(3.21) 

bsTN,<p 

STA CRAT ss 

station clock 

(3.21) 

CSTN,ip 

STA CRRT 88 

model 

(3.21) 


vvvvvv satellite name 
vv satellite number 

as station number 

xm sequence number 
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Table B.I cont. 

Glossary of GPSOMC Parameters 


Parameter 

GPSOMC name 

Definition 

Reference 

f jp 

RSPINAX 

88 

Cylindrical 

(2.1) 

A 

LONGTUD 

88 

station 

(2.2) 

z 

POLPROJ 

88 

coordinates 

(2.3) 

r» p 

DRSP/DT 

88 

Time rates of 

(2.1) 

A 

DLON/DT 

88 

change of 

(2.2) 

z 

DPOL/DT 

88 

station coordinates 

(2.3) 

X 

STA X 

88 

Cartesian 

(2.4) 

Y 

STA Y 

88 

station 

(2.5) 

Z 

STA Z 

88 

coordinates 

(2.6) 

X 

DSTAX/DT 

88 

Time rates of 

(2.4) 

Y 

DSTAY/DT 

88 

change of 

(2.5) 

Z 

DSTAZ/DT 

88 

station coordinates 

(2.6) 

h 

VLOVE 


Vertical Love number 

(2.8) 

l 

HLOVE 


Horizontal Love number 

(2.8) 

*!> 

TIDPHAS 

88 

Tide lag 

(2.13) 

IPP AT 

GEN REL GAMMA 


Post-Newtonian gamma 

(3.8) 

«GC 

GEOCENTER X 


Coordinate frame offset 

(2.79) 

VGC 

GEOCENTER Y 


Cartesian 

(2.79) 

ZGC 

GEOCENTER Z 


components 

(2.79) 

a 

COORD SCALE 


Scale factor 

(2.81) 

ZSVo 

X vvvwv 

Space vehicle epoch 


Vsv o 

Y vvvwv 

position Cartesian 


*SV o 

Z vvvwv 

components 

Sec. 5.4 

*SVo 

DX vvvwv 

Space vehicle epoch 


Vsv 0 

DY vvvwv 

velocity Cartesian 


*SV 0 

DZ vvvwv 

components 



88 station number 

vvvwv satellite name 
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Table B.I cont. 

Glossary of GPSOMC Parameters 


Parameter 

GPSOMC name 

Definition 

Reference 

©i 

X POLE MOTION 

Pole position 

!■ 

©2 

Y POLE MOTION 

components 

■SfBiyi 

UTl - UTC 

UTl MINUS UTC 

UTl - UTC 

(2.30) 

©i 

X POLE RATE 

Pole position 

(2.27) 

©2 

Y POLE RATE 

time rates 

(2.28) 

H 

UTl -UTC RATE 

UTl - UTC time rate 

(2.29) 

*©*,v,*0 

t ROT TWEAK OFFS 

Perturbation rotation 

(2.76) 


t ROT TWEAK RATE 

coefficients 

(2.76) 

Aqj , Alj 

NUT AMPLPSI ynnn 

Nutation amplitudes 

(2.51, 2.53) 

A i j j A Z j 

NUT AMPLPSITynnn 

in longitude 

(2.51, 2.53) 

Aip 

NUT AMPLPSI A 

Longitude nutation tweak 

(2.57) 

Bo j, Bij 

NUT AMPLEPS ynnn 

Nutation amplitudes 

(2.52, 2.54) 

Blj, B 3 y 

NUT AMPLEPS Tynnn 

in obliquity 

(2.52, 2.54) 

Ac 

NUT AMPLEPSA 

Obliquity nutation tweak 

(2.58) 

PZ ir , 

DRYZTROP is 

Dry senith delay 

(4.1) 

PZ„.t 

VETZTROP bb 

Wet senith delay 

(4.1) 

PZir, 

DDTRP/DT as 

Zenith delay 

(4.2) 

Pz m . t 

DWTRP/DT 88 

time rates 

(4.2) 

P 

DRYMAPSG 88 

Lanyi map parameter 

(5.56) 


t X, Y, or Z 

y S or C for sine or cosine 

nnn number of term in Wahr series 

as station number 
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